##################################################################
# File:      BJPS_Appendix_Replication.R
# Purporse:  Replication of Appendix figures and tables for 
#                 British Journal of Political Science
# Author:    Harry Oppenheimer
# Date:      02/27/25
# Last Edit:  
# Data In:   ~/Replication
# Data Out:  ~/Replication/Appendix
# Prev file: None
# Status:    In Progress.
# Computer: 2.6 GHz 6-Core Intel Core i7 laptop with 32 GB of DDR4 RAM, MacOS version 13.5, R 4.2.1
##################################################################

rm(list = ls())

#Load in the libraries for the analysis
libraries <- c(  "data.table","tidyverse", "fixest", "ggpubr","PanelMatch",
                 "igraph", "graphlayouts",
                 "ggraph", "RColorBrewer", "network")
new.packages <- libraries[!(libraries %in% installed.packages()[,"Package"])]
if(length(new.packages)>0){install.packages(new.packages)} 
lapply(libraries, require, character.only = TRUE)

#Set the seed for the analysis
set.seed(02138)

#List the dictionary to make the tables later
table_dictionary = c("allied"= "Treaty",
                     "binary_pta" = "Binary PTA",
                     "wto_joint" = "Joint WTO", 
                     "time_seq" = "Month",
                     "pair" = "Dyad",
                     "paird" = "Directed Dyad",
                     "from" = "CountryA",
                     "to" = "CountryB",
                     "n_agreements_all" = "Interconnection Agreements",
                     "diff_n_agreements_all" = "Differenced Interconnection Agreements",
                     "fromt"= "CountryA * Month",
                     "tot" = "CountryB * Month",
                     "log_agreements_all" = "log(Agreements + 1)")

#Set the working directory to the replication files location
setwd("~/Dropbox/Research/Projects/Replication")
dir.create("Appendix") 

#Load in the data
data <- fread("bjps_data_all.csv")
data_soe <- fread("bjps_data_SOE.csv")

#Log + 1 the DV for the Table 6 analysis 
data$log_agreements_all <- log(data$n_agreements_all + 1)

#Difference the DV for the Table 7 analysis 
data <- data %>%
  group_by(paird) %>%
  arrange(month) %>% 
  mutate(diff_n_agreements_all = n_agreements_all - lag(n_agreements_all))


###############
#Code to create Figure 1 in the Appendix
#Time to run: 1 minute 20 seconds
###############

f <- data %>% group_by(paird) %>%
  summarize(m = mean(allied)) %>%
  filter(m != 0 & m != 1)

data_df <- data %>% group_by(paird) %>%
  mutate(ID = cur_group_id()) %>%
  filter(paird %in% f$paird)

data_df <- as.data.frame(data_df)

data_panel <- PanelData(data_df, unit.id = "ID", time.id = "time_seq",
                        treatment = "allied", outcome = "n_agreements_all")

DisplayTreatment(panel.data = data_panel, legend.position = "none",
                 xlab = "Month", ylab = "Dyad")

ggsave("Appendix/Figure1.png", width = 20, height  = 18,
       device = "png", dpi = 600)

rm(data_df); rm(f); rm(data_panel) 

###############
#Code to create Figure 2 in the Appendix
###############

data %>% group_by(month) %>%
  summarize(prop_allied = n_distinct(pair[allied == 1])/n_distinct(pair)) %>%
ggplot(aes(x = month, y = prop_allied)) + geom_line() +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  ggtitle("Proportion of Dyads With Active ATOP Treaty")

ggsave("Appendix/Figure2.png", width = 8, height  = 4,
       device = "png", dpi = 600)

###############
#Code to create Figure 3 in the Appendix
###############

nwk <- data %>% ungroup %>%
  filter(!duplicated(pair)) %>%
  select(from, to, cable_connect_2010) %>%
  filter(cable_connect_2010 == 1) 

colnames(nwk)[3] <- "weight"
network <- graph.data.frame(nwk, directed = F)

# make degree for labelling most popular nodes
V(network)$degree <- degree(network)

# remove isolates
isolates <- V(network)[degree(network)==0]
network <- delete_vertices(network, isolates)

# calculate modularity for coloring
communities <- cluster_louvain(network)
V(network)$modularity <- communities$membership
members <- membership(communities)

m <- data.frame(table(communities$membership))
m$Var1 <- as.numeric(levels(m$Var1))[m$Var1]
#m <- m[1:which(m$Freq==1)[1],]
m$color <- RColorBrewer::brewer.pal(nrow(m), "Set3")

V(network)$modularity <- ifelse(V(network)$modularity<max(m$Var1), V(network)$modularity, max(m$Var1))
V(network)$color <- m$color[V(network)$modularity]

layout <- create_layout(network, layout = 'fr')
ggraph(layout) +
  geom_node_point(aes(color  =V(network)$color), alpha = .8, size = 10, show.legend = FALSE) +
  geom_edge_link(show.legend = FALSE, alpha = .1)+
  geom_node_text(aes(label = names(V(network))), repel = T, size=3) +
  theme_void()

ggsave("Appendix/Figure3.png", width = 14, height  = 14,
       device = "png", dpi = 600)

rm(m); rm(network); rm(nwk); rm(layout); rm(communities)

###############
#Code to create Table 6 in the Appendix
###############

m1 = feols(log_agreements_all ~ allied | paird  +  fromt  + tot,
                 data,
                 cluster = ~paird+time_seq)

m2 = feols(log_agreements_all ~ allied  + wto_joint | paird  +  fromt  + tot,
                     data,
                     cluster = ~paird+time_seq)

m3 = feols(log_agreements_all ~ allied  + binary_pta | paird  +  fromt  + tot,
                       data,
                       cluster = ~paird+time_seq)

m4 = feols(log_agreements_all ~ allied  + binary_pta + wto_joint | paird  +  fromt  + tot,
                     data,
                     cluster = ~paird+time_seq)

etable(m1, m2,  m3, m4, tex = T, digits = 4, dict = table_dictionary)


###############
#Code to create Table 7 in the Appendix
###############

m1 = feols(diff_n_agreements_all ~ allied | fromt  + tot,
                 data,
                 cluster = ~paird+time_seq)

m2 = feols(diff_n_agreements_all ~ allied  + wto_joint |  fromt  + tot,
                     data,
                     cluster = ~paird+time_seq)

m3 = feols(diff_n_agreements_all ~ allied  + binary_pta | fromt  + tot,
                       data,
                       cluster = ~paird+time_seq)

m4 = feols(diff_n_agreements_all ~ allied  + binary_pta + wto_joint | fromt  + tot,
                     data,
                     cluster = ~paird+time_seq)

etable(m1, m2,  m3, m4, tex = T, digits = 4, dict = table_dictionary)

rm(m1); rm(m2); rm(m3); rm(m4)

###############
#Code to create Figure 4 in the Appendix
###############

f <- data %>% group_by(pair) %>%
  summarize(m = mean(allied)) %>%
  filter(m != 0 & m != 1)

treated_pairs <- as.vector(f$pair)

coef_df <- data.frame(Pair = treated_pairs, Estimate = NA, se = NA)

for(i in 1:length(treated_pairs)){
  datatemp <- data[-which(data$pair == treated_pairs[i]),]
  
  m1 = feglm(n_agreements_all ~ allied | fromt + tot  + paird,
              datatemp, family=quasipoisson, notes = FALSE)
  
  coef_df[i,2] <- m1$coeftable[1]
  coef_df[i,3] <- se(m1)
}


coef_df %>%
  mutate(a = ifelse(Estimate - 1.96*se < 0, 0, 1)) %>%
  ggplot() + geom_point(aes(x = Pair, y = Estimate)) + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + 
  geom_linerange(aes(x = Pair, 
                     ymin = Estimate - 1.96*se,
                     ymax = Estimate + 1.96*se,
                     colour = factor(a)),
                 lwd = .5) + 
  scale_x_discrete(guide = guide_axis(n.dodge=3)) +
  theme(legend.position = "none", axis.text.x = element_text(angle = 45)) + 
  xlab("Dyad")

ggsave("Appendix/Figure4.png", width = 12, height  = 8,
       device = "png", dpi = 600)


